#include <time.h>
#ifndef SERIALQGCM
#include <omp.h>
#endif

#include "mappers.h"
#include "qgcmmath.h"
/* ------------------------------------------------------------------------- */
custom_func_container::custom_func_container(const string & name, const int nangs, const int nks, const int nfs)
	: 
		mo(getModel()),
		OPname(name),
		rank_fields(nfs),
		rank_kspace(nks),
		rank_amodes(nangs)
{
	for(int i = 0; i < rank_fields; i++){
		Value * pval = new Value( static_cast<int_type>(0) );
		fields.push_back( pair<string, Value *>( "", pval) );
	}
	for(int i = 0; i < rank_kspace; i++){
		Value * pval = new Value( static_cast<float_type>(0) );
		kvariables.push_back( pair<string, Value *>( "", pval) );
	}
	for(int i = 0; i < rank_kspace; i++){
		Value * pval = new Value( static_cast<int_type>(0) );
		angmodevariables.push_back( pair<string, Value *>( "", pval) );
	}
	cval = NULL;
	/*cval.DefineConst( "R_coeff" , static_cast<float_type>(-3117.09091) ); // factor of -2*(2*pi)^4
	  cval.DefineConst( "S_coeff" , static_cast<float_type>(-6234.18183) ); // factor of -4*(2*pi)^4*/
}
/* ------------------------------------------------------------------------- */
void custom_func_container::assign_angspace_names( const vector<string> & names )
{
	rassert( names.size() == rank_amodes );
	for(int i = 0; i < rank_amodes; i++){
		angmodevariables[i].first = names[i];
	}
}
/* ------------------------------------------------------------------------- */
void custom_func_container::assignkspacenames( const vector<string> & names ){
	rassert( names.size() == rank_kspace );
	for(int i = 0; i < rank_kspace; i++){
		kvariables[i].first = names[i];
	}
}
/* ------------------------------------------------------------------------- */
void custom_func_container::assign_field_names( const vector<string> & names ){
	rassert( names.size() == rank_fields );
	for(int i = 0; i < rank_fields; i++){
		fields[i].first = names[i];
	}
}
/* ------------------------------------------------------------------------- */
void custom_func_container::assignindices(const rank2 & r1 ){
	for( int i = 0; i < 2 && i < this->fields.size(); i++){
		*(this->fields[i].second) = static_cast<int_type>( r1[i] );
	}
}
/* ------------------------------------------------------------------------- */
void custom_func_container::assignindices(const rank3 & r1 ){
	for( int i = 0; i < 3 && i < this->fields.size(); i++){
		*(this->fields[i].second) = static_cast<int_type>( r1[i] );
	}
}
/* ------------------------------------------------------------------------- */
int custom_func_container::getnum_ang_modes() const {
	return this->rank_amodes;
};
/* ------------------------------------------------------------------------- */
int custom_func_container::getnum_kspace_variables() const {
	return this->rank_kspace;
};
/* ------------------------------------------------------------------------- */
int custom_func_container::getnum_fields() const {
	return this->rank_fields;
};
/* ------------------------------------------------------------------------- */
void custom_func_container::assignindices(const rank3 & r1, const rank3 & r2 ){
	for( int i = 0; i < 3 && i < this->fields.size(); i++){
		*(this->fields[i].second) = static_cast<int_type>( r1[i] );
	}
	for( int i = 3; i < 6 && i < this->fields.size(); i++){
		*(this->fields[i].second) = static_cast<int_type>( r2[i-3] );
	}
}
/* ------------------------------------------------------------------------- */
void custom_func_container::assignnkspace( INDEX nk, INDEX np ){
        assert( this->getnum_ang_modes() == 2 );
	*(this->angmodevariables[0].second) = static_cast<int_type>(nk);
	*(this->angmodevariables[1].second) = static_cast<int_type>(np);
}
/* ------------------------------------------------------------------------- */
void custom_func_container::assignkspace(REAL k1){
        assert( this->getnum_kspace_variables() == 1 );
	*(this->kvariables[0].second) = static_cast<float_type>(k1);
}
/* ------------------------------------------------------------------------- */
void custom_func_container::assignkspace(REAL k1, REAL k2){
        assert( this->getnum_kspace_variables() == 2 );
	*(this->kvariables[0].second) = static_cast<float_type>(k1);
	*(this->kvariables[1].second) = static_cast<float_type>(k2);
}
/* ------------------------------------------------------------------------- */
void custom_func_container::assignkspace(REAL k1, REAL k2, REAL k3){
        assert( this->getnum_kspace_variables() == 3 );
	*(this->kvariables[0].second) = static_cast<float_type>(k1);
	*(this->kvariables[1].second) = static_cast<float_type>(k2);
	*(this->kvariables[2].second) = static_cast<float_type>(k3);
}
/* ------------------------------------------------------------------------- */
const IValue & custom_func_container::eval(){
	return this->cval->Eval();
}
/* ------------------------------------------------------------------------- */
const string & custom_func_container::name() const{
	return this->OPname;
}
/* ------------------------------------------------------------------------- */
custom_func_container::~custom_func_container(){
	delete cval;
	for(int i = 0; i < fields.size(); i++) delete fields[i].second;
	for(int i = 0; i < kvariables.size(); i++) delete kvariables[i].second;
	for(int i = 0; i < angmodevariables.size(); i++) delete angmodevariables[i].second;
}
/* ------------------------------------------------------------------------- */
bool custom_func_container::isdefined() const {
	return this->cval != NULL;
}
/* ------------------------------------------------------------------------- */
void custom_func_container::assignexpr(const token_list & expr){
	// failure to evaluate expression may lead to an exception being
	// thrown.
	delete [] this->cval;
	// allocate the parser
	this->cval = new ParserX(pckNONE); // testing parser
	mo.assignparserdata(*cval, "ifc", false);
	mo.assigndiagfunctions(*cval);
	// assign all variables
	for(int i = 0; i < rank_fields; i++){
		rassert( fields[i].first.size() > 0 );
		//cval->DefineVar( fields[i].first, Variable(* (fields[i].second) ) );
		forceDefineVar( *cval, fields[i].first, fields[i].second );
	}
	for(int i = 0; i < rank_kspace; i++){
		rassert( kvariables[i].first.size() > 0 );
		//cval->DefineVar( kvariables[i].first, Variable(* (kvariables[i].second) ) );
		forceDefineVar( *cval, kvariables[i].first, kvariables[i].second );
	}
	for(int i = 0; i < rank_amodes; i++){
		rassert( angmodevariables[i].first.size() > 0 );
		//cval->DefineVar( angmodevariables[i].first, Variable(* (angmodevariables[i].second) ) );
		forceDefineVar( *cval, angmodevariables[i].first, angmodevariables[i].second );
	}
	this->cval->SetExpr(expr);
	// try to parser the expression 
	this->cval->Eval(); // exception should be handled by caller
}
// =================================================================================================
custom_func_container get_custom_func_container::asym_kpq_cp("asymkpq", 0, 3, 0);
custom_func_container get_custom_func_container::sym_kpq_cp("symkpq", 2, 3, 0);
custom_func_container get_custom_func_container::R_cp("R", 2, 3, 6);
custom_func_container get_custom_func_container::S_cp("S", 2, 3, 6);
custom_func_container get_custom_func_container::customAhat("A", 0, 3, 3);
custom_func_container get_custom_func_container::customQ("Q", 0, 1, 2);
custom_func_container get_custom_func_container::customP("P", 0, 0, 3);


get_custom_func_container::get_custom_func_container(const string & name){
	wfunc = NULL;
	if(asym_kpq_cp.name() == name){
		wfunc =  &(asym_kpq_cp);
	}else if(R_cp.name() == name){
		wfunc =  &(R_cp);
	}else if(S_cp.name() == name){
		wfunc =  &(S_cp);
	}else if(sym_kpq_cp.name() == name){
		wfunc =  &(sym_kpq_cp);
	}else if(customAhat.name() == name){
		wfunc =  &(customAhat);
	}else if(customQ.name() == name){
		wfunc =  &(customQ);
	}else if(customP.name() == name){
		wfunc =  &(customP);
	}
}

get_custom_func_container::get_custom_func_container(int index){
	wfunc = NULL;
	if(index == 1){
		wfunc =  &(asym_kpq_cp);
	}else if(index == 2){
		wfunc =  &(R_cp);
	}else if(index == 3){
		wfunc =  &(S_cp);
	}else if(index == 4){
		wfunc =  &(sym_kpq_cp);
	}else if(index == 5){
		wfunc =  &(customAhat);
	}else if(index == 6){
		wfunc =  &(customQ);
	}else if(index == 7){
		wfunc =  &(customP);
	}else{	
		fatal_stop("Bad argument provided to container function.");
	}
}

custom_func_container & get_custom_func_container::operator*() {
	return * wfunc;
}

bool get_custom_func_container::isvalid() const {
	return wfunc != NULL;
}

bool get_custom_func_container::operator!() const {
	return wfunc == NULL;
}
// =================================================================================================
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void custom_symcomputeKPQintegrand(REAL & rvalue, const REAL & k, const REAL & p, const REAL & q){
	custom_func_container & cp = *get_custom_func_container(4);
	cp.assignkspace(k,p,q);
	rvalue = static_cast<REAL>( cp.eval().GetFloat() );
};
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
qgcomplex custom_asymcomputeKPQintegrand(const REAL & k, const REAL & p, const REAL & q, const INDEX & np, const INDEX & mk)
{
	custom_func_container & cp = *get_custom_func_container(1);  
	cp.assignkspace(k,p,q);
	cp.assignnkspace(mk, np);
	const cmplx_type & cval = cp.eval().GetComplex();
	return qgcomplex( cval.real(), cval.imag() );
	//return qgcomplex( static_cast<REAL>(cp.freal.Eval()), static_cast<REAL>(cp.fimag.Eval() ) );
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
REAL Rcustom(const model_class & mo, REAL k2, REAL p2, REAL q2, size_t i1, size_t i2){
	custom_func_container & cp = *get_custom_func_container(2);
	double k = sqrt(k2); 
	double p = sqrt(p2);
	double q = sqrt(q2);
	cp.assignkspace( k, p, q );
	const vector<rank3> & aindices = mo.getAfieldindices();
	cp.assignindices( aindices[i1], aindices[i2] );
	return static_cast<REAL>( cp.eval().GetComplex().real() );
};
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
REAL Scustom(const model_class & mo, REAL k2, REAL p2, REAL q2, size_t i1, size_t i2){
	custom_func_container & cp = *get_custom_func_container(3);
	double k = sqrt(k2); 
	double p = sqrt(p2);
	double q = sqrt(q2);
	cp.assignkspace( k, p, q );
	const vector<rank3> & aindices = mo.getAfieldindices();
	cp.assignindices( aindices[i1], aindices[i2] );
	return static_cast<REAL>( cp.eval().GetComplex().real() );
};
//. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
REAL customQ(const model_class &, const REAL & k, const rank2 & field_indices){
	custom_func_container & cp = *get_custom_func_container(6);
	cp.assignkspace(k);
	rank3 rpos( field_indices[0], field_indices[1], field_indices[1] );
	cp.assignindices( rpos, rpos );
	return static_cast<REAL>( cp.eval().GetComplex().real() );
};
//. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
REAL customAhatfunc(const model_class & mo, const REAL & k, const REAL & p, const REAL & q, const rank3 & ipos){
	custom_func_container & cp = *get_custom_func_container(5);
	cp.assignkspace( k, p, q);
	cp.assignindices( ipos, ipos ); 
	return static_cast<REAL>( cp.eval().GetComplex().real() );
}; 
//. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
bool customP(const rank3 & ipos){
	custom_func_container & cp = *(get_custom_func_container(7));
	cp.assignindices( ipos, ipos );
	if(cp.eval().GetComplex().real() != 0) return true;
	return false;
}; 

// --------------------------------------------------------------- //
// class mirror_array definition:
mirror_array::mirror_array(const vector<INDEX> & array){
	this->data = new vector<INDEX>;
	*(this->data) = array;
};

const INDEX & mirror_array::back() const{
	if(this->data){
		return this->data->back();
	}else{
		fatal_stop("No array was assigned to container class.");
	}
	return this->data->front();
};

mirror_array::mirror_array(){
	this->data = new vector<INDEX>;
};
mirror_array::~mirror_array(){
	delete this->data;
	this->data = NULL;
};

mirror_array & mirror_array::operator=(const vector<int> & other ){
	vector<INDEX> & array = *(this->data);
	array.clear();
	array.resize(other.size(), 0);
	for(size_t i = 0; i < other.size(); i++){
		array[i] = static_cast<INDEX>(other[i]);
	}
	return *this;
};
size_t mirror_array::size() const{
	if(this->data){
		return this->data->size();
	}
	fatal_stop("No array was assigned to container class.");
	return 0;
}
INDEX  mirror_array::operator[](INDEX index) const{
	if(this->data){
		if(index < 0){
			size_t i = static_cast<size_t>( -1*index );
			return -1*( (*data)[i]);
		}else{
			return (*data) [ static_cast<size_t> ( index ) ] ;
		}
	}else{
		fatal_stop("No array was assigned to container class.");
		return 0;
	}
};
//----------------------------------------------------------------------------//
// This function is much faster in fortran, and is only used for debugging here
REAL computeLawSinesCoeff(REAL k, REAL p, REAL q){
	REAL k2 = k*k; REAL p2 = p*p; REAL q2 = q*q;
	return (1.0/(2.0*k*p*q)) * sqrt( (k2 + p2 + q2)*(k2+p2+q2) -2.0*(k2*k2 + p2*p2 + q2*q2));
}
/* ----------------------------------------------------------------------------*/
inline bool computestaticpart::isinsideregularizedslot(const REAL k, const REAL p, const REAL q) const{
	if(ignoreslot) return true;
	if(ignorebasetriangle){
		if( q < kbasefactor*k && p < kbasefactor*k ) return false;
	}
	return ( q < k && p <= k + q && p >= k - q ) || ( q >= k && p <= k + q && p >= -k + q);
};
/* ----------------------------------------------------------------------------*/
inline bool computestaticpart::isinsideslot(const REAL k, const REAL p, const REAL q) const{
	/*if(ignorebasetriangle){
	  if( q < kbasefactor*k && p < kbasefactor*k ) return false;
	  }*/
	if(ignoreslot) return true;
	return ( q < k && p <= k + q && p >= k - q ) || ( q >= k && p <= k + q && p >= -k + q);
}; 
/* ----------------------------------------------------------------------------*/
inline REAL * calculate_midpoints(const REAL xmin, const REAL dx, const int nresol, bool log_spacing){
	REAL * midpoints = new REAL[nresol];
	if(log_spacing){
		REAL xi = xmin;
		for(int i = 0; i < nresol; i++){
			midpoints[i] = dx/(log(1.0+dx/xi));
			xi+=dx;
		}
	}else{
		REAL xf = xmin + 0.5*dx;
		for(int i = 0; i < nresol; i++){
			midpoints[i] = xf;
			xf += dx;
		}
	}
	return midpoints;
}
/* ----------------------------------------------------------------------------*/
REAL computestaticpart::fine_mesh_integrate(const INDEX ki, const INDEX pi, const INDEX qi,
		const float * width_factors) const
{
	const model_class & model = *mo;
	const bool k_log_spacing = kspacing[ki] != gdb_enum_UNIFORMSPACING;
	const bool p_log_spacing = kspacing[pi] != gdb_enum_UNIFORMSPACING;
	const bool q_log_spacing = kspacing[qi] != gdb_enum_UNIFORMSPACING;
	const REAL kc = kdomain_centered[ki];
	const REAL pc = kdomain_centered[pi];
	const REAL qc = kdomain_centered[qi];
	{
		model.printfmsg("Fine mesh integrating (k,p,q)=( %f, %f, %f)\n", kc,pc,qc);
		model.flushmessages();
	}
	REAL lvol = 1.0; // last computed volume
	REAL vol = 0.0;
	const size_t ksize = kdomain_centered.size();
	const REAL kmax = kdomain[ki+1];
	const REAL kmin = kdomain[ki];
	const REAL pmax = kdomain[pi+1];
	const REAL pmin = kdomain[pi];
	const REAL qmax = kdomain[qi+1];
	const REAL qmin = kdomain[qi];
	const REAL widthp = pmax - pmin;
	const REAL widthq = qmax - qmin;
	const REAL widthk = kmax - kmin;

	bool firstloop = true;
	unsigned int depth = 0; // nesting depth, 0 refers to the first fine domain
	// in memory cases depth = 0 may not be available
	int npresol, nqresol, nkresol;
	// we use as a first guess func evaluated at the center * volume of the box
	symkpqintegrand(lvol, kc, pc, qc); // midpoint evaluation of the integrand 
	lvol*= (widthp*widthq*widthk);
	REAL err_real = 1; // = fabs(vol-lvol)/fabs(lvol);
	while( firstloop || err_real > THRESHOLD )
	{ // fine mesh integration loop
		if(depth < nresols.size()){ 
			int nresol = nresols[depth];
			unsigned int kpos = (static_cast<int>(ki)*ksize+static_cast<int>(pi))*ksize + static_cast<int>(qi);
			unsigned int ppos = (static_cast<int>(pi)*ksize+static_cast<int>(qi))*ksize + static_cast<int>(ki);
			unsigned int qpos = (static_cast<int>(qi)*ksize+static_cast<int>(ki))*ksize + static_cast<int>(pi);
			assert(kpos < ksize*ksize*ksize);
			assert(ppos < ksize*ksize*ksize);
			assert(qpos < ksize*ksize*ksize);
			const float kfactor = width_factors[kpos];
			const float pfactor = width_factors[ppos];
			const float qfactor = width_factors[qpos];
			nkresol = static_cast<int>(static_cast<float>(nresol)*kfactor);
			if(nkresol < 1) nkresol = 1;
			npresol = static_cast<int>(static_cast<float>(nresol)*pfactor);
			if(npresol < 1) npresol = 1;
			nqresol = static_cast<int>(static_cast<float>(nresol)*qfactor);
			if(nqresol < 1) nqresol = 1;
		}else{
			// we have a problem not enough deptbh for the system to converge
			// add parallel support here ...
			model.printfmsg("Symmetric static integration did not converge under finest mesh resolution.\n");
			model.printfmsg("Please increase the size of the multiplier array.\n");
			model.flushmessages();
			throw cgcmexception("Error in computation of interaction coefficients.", -1);
		}
		model.printfmsg("new mesh, depth='%i' sampling (numk,nump,numq) = (%i,%i,%i)",
				depth, nkresol, npresol, nqresol);
		model.printfmsg(".\n");
		if(!firstloop){
			lvol = vol; // for convergence tests
			vol = 0;
		}
		// we increment pf**2 instead of pf and qf**2 instead of qf
		// to save calculations of the squares of the k-space vector magnitudes here
		REAL dkf, dpf, dqf;
		REAL * kf_midpoints, * pf_midpoints, * qf_midpoints;
		{ 
			dkf=widthk/static_cast<REAL>(nkresol); 
			dpf=widthp/static_cast<REAL>(npresol); 
			dqf=widthq/static_cast<REAL>(nqresol);
			// first fill in the mid points for k
			kf_midpoints = calculate_midpoints( kmin, dkf, nkresol, k_log_spacing );
			if(pi == ki) pf_midpoints = kf_midpoints;
			else  pf_midpoints = calculate_midpoints( pmin, dpf, npresol, p_log_spacing );
			if(qi == ki) qf_midpoints = kf_midpoints;
			else if(qi == pi) qf_midpoints = pf_midpoints;
			else qf_midpoints = calculate_midpoints( qmin, dqf, nqresol, q_log_spacing );
			for(int kcounter = 0; kcounter < nkresol; kcounter++){
				const REAL & kf = kf_midpoints[kcounter];
				for(int pcounter = 0; pcounter < npresol; pcounter++){
					const REAL & pf = pf_midpoints[pcounter];
					for(int qcounter = 0; qcounter < nqresol; qcounter++){
						const REAL & qf = qf_midpoints[qcounter];
						if( isinsideslot(kf,pf,qf) ){ //makfe sure we are inside the slot domain ....
							REAL fvol; symkpqintegrand(fvol, kf, pf, qf);
							if(fvol == fvol) vol += fvol;
						}
					}}}
					vol*= (dkf*dpf*dqf);
					firstloop = false;
					depth++;
					err_real = fabs(vol-lvol)/fabs(lvol);
					model.printfmsg("---> mesh completed with new val=%f, last val=%f, rel error=%f.\n", vol, lvol, err_real);
					delete [] kf_midpoints;
					if(ki != pi ) delete [] pf_midpoints;
					if(qi != ki && qi != pi) delete [] qf_midpoints;
		} // close fine mesh integration loop
	}
	model.printfmsg( "Converged @ depth, error = %i, %f.\n", depth, err_real);
	model.printfmsg( "- - - - - - - - - - - - - - - - - - - - - - -\n"); 
	model.printfmsg( "Fine mesh converged with value = %f.\n", vol);
	model.printfmsg( "---------------------------------------------\n"); 
	model.flushmessages();
	return vol;
}
/* ----------------------------------:-------------------------------------------------------------------------------------*/
rank3 computestaticpart::sortrank3(const rank3 & pos){
	// this function orders the indices from lowest to highest
	INDEX a = pos[0];
	INDEX b = pos[1];
	INDEX c = pos[2];
#define swap(x, y) {INDEX tmp; tmp = x; x = y; y = tmp; } 
	if (b < a) swap(b, a);
	if (c < b) swap(c, b);
	if (b < a) swap(b, a);
#undef swap
	return rank3(a,b,c);
};
/* ------------------------------------------------------------------------------------------------------------------------*/
// This function goes through and assigns either 0 if the value is to be calculated later using
// fine_mesh_integrate or
// the value of the index in the array this specific iter should be assigned to if the value
// has been previously calculated in the array
int computestaticpart::findsymvalue(const rank3 * symmetric_data_map, const rank3 & pos, const size_t & miter)const{
	rank3 sorted_pos = sortrank3(pos); 
	for(size_t iter = 0; iter < miter; ++iter){
		const rank3 & cpos = symmetric_data_map[iter];
		if( cpos == sorted_pos ){
			return iter;
		}
	}
	return -1;
};

/* ------------------------------------------------------------------------------------------------------------------------*/
bool computestaticpart::preassignsymvalue(
		REAL & sympart, const REAL * symmetric_data, const rank3 * symmetric_data_map, const rank3 & pos, const size_t & miter)
const {
	// first sort the index that was given
	rank3 sorted_pos = sortrank3(pos); 
	for(size_t iter = 0; iter < miter; ++iter){
		const rank3 & cpos = symmetric_data_map[iter];
		if( cpos == sorted_pos ){
			sympart = symmetric_data[iter]; 
			return true;
		}
	}
	return false;
};
#ifdef DEBUG
/* ------------------------------------------------------------------------------------------------------------------------*/
void dumpsymmetricdata(REAL * data, rank3 * datamap, int size){
	FILE * fout = getrout();
	for(int i = 0; i < size; i++){
		const rank3 & r3 = datamap[i];
		fprintf(fout, "data[%i,%i,%i] = %g\n", static_cast<int>(r3[0]),
				static_cast<int>(r3[1]),
				static_cast<int>(r3[2]), data[i]);
	} 
};
#endif
/* --------------------------------------------------------------------------------------------------------------*/
computestaticpart::computestaticpart(
		const model_class & model, 
		static_container * & Rs,
		static_container * & Ss,
		map_container<REAL, 6> * & Amap,
		map_container<REAL, 3> * & Qmap,
		map_container<REAL, 3> * & QImap,
		void ( * kpqintegrand_arg ) ( REAL &, const REAL &, const REAL &, const REAL & ),
		complex <REAL> ( * asymkpqintegrand_arg )( const REAL & , const REAL &, const REAL &, const INDEX &, const INDEX &),
		bool ( * fcompute_P ) ( const rank3 & ),
		REAL ( * fcompute_Q ) ( const model_class & mo, const REAL &, const rank2 & ),
		REAL ( * fcompute_A ) ( const model_class & mo, const REAL &, const REAL &, const REAL &, const rank3 &),
		const REAL R_coeff_arg,
		const REAL S_coeff_arg
		)   : R_coeff(0), S_coeff(0), ignoreslot(false), NUMPROCESSORS(1), numfields(1), THRESHOLD(0.01), ignorebasetriangle(true),  dk_width_adjuster(0.25), kbasefactor(1.0)
{
	this->mo = &(model);
	this->Rstaticweights = new static_container;
	this->Sstaticweights = new static_container;
	this->Ahatmap = new map_container<REAL, 6>;
	this->Qmap    = new map_container<REAL, 3>;
	this->QImap   = new map_container<REAL, 3>;
	// 4 function pointers required for the calculation
	this->func_compute_P = fcompute_P;
	this->func_compute_Q = fcompute_Q;
	this->func_compute_A = fcompute_A;
	this->asymkpqintegrand = asymkpqintegrand_arg;
	this->symkpqintegrand = kpqintegrand_arg;
	*( const_cast<REAL *>(&R_coeff)) = R_coeff_arg;
	*( const_cast<REAL *>(&S_coeff)) = S_coeff_arg;
	this->initialize();
	this->mapstaticpart();
	Rs = this->Rstaticweights;
	Ss = this->Sstaticweights;
	Amap = this->Ahatmap;
	Qmap = this->Qmap;
	QImap = this->QImap;
};
#ifdef DEBUG
/* --------------------------------------------------------------------------------------------------------------*/
void computestaticpart::initialize_symmetric_data_array(const size_t symmetric_data_size, REAL *symmetric_data,
		rank3 * symmetric_data_map ){
	{ 
		const REAL fNaN = sqrt(-1.0);
		for(size_t iter = 0; iter < symmetric_data_size; iter++){
			symmetric_data[iter] = fNaN;  
			symmetric_data_map[iter] = rank3(-1,-1,-1);
		}
	}
};
#endif
/* --------------------------------------------------------------------------------------------------------------*/
void computestaticpart::initialize(){
#define assign_param(type,arg1,arg2) *(const_cast<type*>(&(this->arg1)))=arg2 ;
	const model_class & model = * mo;
	assign_param(REAL,THRESHOLD,model.get_global_param_by_name(gdb_str_PQFINEDOMAINACCURACYTHRESHOLD).getfValue());
	assign_param(REAL,dk_width_adjuster,model.get_global_param_by_name(gdb_str_SLOT_WIDTH_SCALE).getfValue());
	assign_param(bool,ignorebasetriangle,model.get_global_param_by_name(gdb_str_IGNORETRIANGLEBASE).getiValue() != 0);
	assign_param(REAL, kbasefactor, model.get_global_param_by_name(gdb_str_KCUTOFFMULTIPLIER).getfValue()); 
	// . . . . . . . . . . . . . . . . .
	// Retrieve model parameters needed to define the characteristics of the fine domain integration.  
#ifndef SERIALQGCM
	assign_param(int,NUMPROCESSORS,model.get_num_threads());
#endif
	// Load the number of field indices
	assign_param(int,numfields,model.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue());
	assign_param(bool,ignoreslot,model.get_global_param_by_name(gdb_str_IGNORESLOT).getiValue() != 0);
	// Load constant profiles from model
	assign_param(vector<REAL>,kdomain,model.get_global_prof_by_name(gdb_prof_str_K_DOMAIN).getData());
	assign_param(vector<REAL>,dkdomain,model.get_global_prof_by_name(gdb_prof_str_DK_DOMAIN).getData());
	assign_param(vector<REAL>, kdomain_centered, model.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData());
	assign_param(vector<REAL>, kspacing, model.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_SPACING).getData());
	const vector<int> & nkdom = model.get_global_prof_by_name(gdb_prof_str_NK_DOMAIN).getiData();
	*const_cast<mirror_array*>(&(this->nkdomain)) = nkdom;
	const vector<REAL>& vnresols= model.get_global_prof_by_name(gdb_prof_str_FINEKDOMAINMULTIPLIERS).getData();
	vector<int> & resols = *(const_cast< vector<int> * >(&(this->nresols)) );
	resols.resize(vnresols.size(), 0);
	for(size_t i = 0; i < vnresols.size(); i++){
		resols[i] = static_cast<int>( vnresols[i] );
	};

	REAL fourier_factor = model.get_global_param_by_name(gdb_str_FOURIER_TRANSFORM_COEFFICIENT).getfValue();
	set_fourier_transform_coefficient( fourier_factor );

#ifndef SERIALQGCM
	omp_set_num_threads(NUMPROCESSORS);
#endif
}; 
/* --------------------------------------------------------------------------------------------------------------*/
void computestaticpart::mapsympart(const size_t symmetric_data_size, REAL * symmetric_data, rank3 * symmetric_data_map)
	const
{
	// we use the map below to determine which elements 
	//   -1 : to calculate
	//   >0 : to assign using the value within
	//  <-1 : to ignore (NaN)
	int * pre_map = new (nothrow) int[symmetric_data_size];
	rassert(pre_map != NULL);
	for(size_t i = 0; i < symmetric_data_size; i++) pre_map[i] = -2;
#ifdef DEBUG
	this-> initialize_symmetric_data_array( symmetric_data_size, symmetric_data, symmetric_data_map );
#endif
	size_t iter = 0;
	const unsigned int ksize = this->kdomain_centered.size(); 
	float * width_factors = new (nothrow) float[ksize*ksize*ksize];
	rassert( width_factors != NULL );
	// map the width_factors first
	{
		for(INDEX ki = 0; ki < static_cast<INDEX>(ksize); ki++){
			const REAL & dk = dkdomain[ki];
			for(INDEX pi = 0; pi < static_cast<INDEX>(ksize); pi++){ // loop over p
				const REAL & dp = dkdomain[pi];
				for(INDEX qi = 0; qi < static_cast<INDEX>(ksize); qi++){ // loop over q
					const REAL & dq = dkdomain[qi];
					unsigned int pos = (static_cast<int>(ki)*ksize+static_cast<int>(pi))*ksize + static_cast<int>(qi);
					width_factors[pos] = pow(((dk*dk)/(dp*dq)), this->dk_width_adjuster);
					assert(pos < ksize*ksize*ksize);
				}}}
	}
	// note that the searchmap function expects the k,p,q indexmap to be sorted when initially added.
	// that is 0,0,0 then 0,0,1, then 0,1,1 ... 1,0,0, 1,0,1 etc ...
	int num_to_compute = 0;
	for(INDEX ki = 0; ki < static_cast<INDEX>(ksize); ki++){  
		const REAL & k = kdomain_centered[ki];
		for(INDEX pi = 0; pi < static_cast<INDEX>(ksize); pi++){ // loop over p
			const REAL & p = kdomain_centered[pi];
			for(INDEX qi = 0; qi < static_cast<INDEX>(ksize); qi++){ // loop over q
				const REAL & q = kdomain_centered[qi];
				const rank3 pos(ki,pi,qi);
				if( isinsideslot(k,p,q) )
					// only calculate elements whose center lies within the slot
				{
					// this init function assigns the value -1 if the function fine_mesh_integrate is to be
					// called or the value of the array iterator where the value of the current element
					// is to be found.
					int piter = this->findsymvalue(symmetric_data_map, pos, iter);
					assert( iter < symmetric_data_size);
					if(piter >= 0){
						pre_map[iter] = piter; // simply assiging piter works but is less clear
					}else{
						pre_map[iter] = -1; // we need to calculate it using the fine_mesh_integrate
						num_to_compute++;
					}
					symmetric_data_map[iter] = pos;
				}
				iter++;
			}}
	}
	int * iter_map = new (nothrow) int[num_to_compute]; // holds indices in symmetric data map of elements we need to compute
	rassert(iter_map != NULL);
	size_t i_iter = 0;
	for(iter = 0; iter < symmetric_data_size; iter++){
		if(pre_map[iter] == -1){
			iter_map[i_iter] = static_cast<int>(iter);
			assert(i_iter < static_cast<size_t>(num_to_compute));
			i_iter++;
		}
	}
	// had to split the above into separate looops for better parallelization efficiency
#ifndef SERIALQGCM
#pragma omp parallel
#pragma omp for schedule(dynamic, 1)
#endif
	for(i_iter = 0; i_iter < static_cast<size_t>(num_to_compute); i_iter++){
		const int iter = iter_map[i_iter];
		const rank3 & pos = symmetric_data_map[iter];
		const INDEX ki = pos[0];
		const INDEX pi = pos[1];
		const INDEX qi = pos[2];
		REAL fval = this->fine_mesh_integrate(ki, pi, qi, width_factors);
		symmetric_data[iter] = fval;
		assert(iter < static_cast<size_t>(symmetric_data_size));
	}
	delete [] iter_map;
	delete [] width_factors;
	mo->printfmsg("Assigning all symmetric values to symmetric data set.\n");
	// now go through and assign all assignment values
	for(iter = 0; iter < symmetric_data_size; iter++){
		if(pre_map[iter] >= 0){
			symmetric_data[iter] = symmetric_data[pre_map[iter]];
		}
	}
	delete [] pre_map;
};
/* --------------------------------------------------------------------------------------------------------------*/
inline int getAindexfromset(
		const unsigned int & ksize_arg, const int & numfields, const INDEX ki_arg, const INDEX pi_arg, const INDEX qi_arg,
		const INDEX ii, const INDEX jj, const INDEX kk
		){
	int ksize  = static_cast<int>(ksize_arg);
	int iindex = static_cast<int>(ii);
	int jindex = static_cast<int>(jj);
	int kindex = static_cast<int>(kk);
	int ki = static_cast<int>(ki_arg);
	int pi = static_cast<int>(pi_arg);
	int qi = static_cast<int>(qi_arg);
	int findex = kindex + numfields*(jindex + numfields*(iindex + numfields*(qi+ksize*(pi+ksize*ki))));
	return findex;
}; 
/* --------------------------------------------------------------------------------------------------------------*/
int findAiter(const vector<rank3> & list, const rank3 & pos){
	for(size_t i = 0; i < list.size(); i++){
		if( list[i] == pos ) return static_cast<int>(i);
	}
	return -1;
};
/* --------------------------------------------------------------------------------------------------------------*/
#ifdef DEBUG
void multiplymatrix(vector< vector<double> > & prod, const vector< vector<double> > & m1, const vector< vector<double> > & m2){
	assert( m1.size() == m2.size() );
	assert( m2.size() > 0 );
	assert( m1.size() == m2[0].size() );
	assert( m1[0].size() == m2.size() );
	prod.resize( m1.size(), vector<double>( m2[0].size(), 0) ); 
	for(size_t i = 0; i < m1.size(); i++){
		for(size_t j = 0; j < m2[0].size(); j++){
			float val = 0;
			for( size_t m = 0; m < m1[0].size(); m++){
				val += m1[i][m]*m2[m][j];
			}
			prod[i][j] = val;
		}}
};
/* --------------------------------------------------------------------------------------------------------------*/
void dumpmatrix(const vector< vector<double> > & matrix){
	FILE * rout = getrout();
	fprintf(rout, "{\n");
	for(size_t i = 0; i < matrix.size(); i++){
		for(size_t j = 0; j < matrix[i].size(); j++){
			if(j > 0){
				fprintf(rout, ", ");
			}else{ 
				fprintf(rout, "{ ");
			}
			fprintf(rout, "%15g", matrix[i][j]);
		}
		fprintf(rout, " }");
		if( i < matrix.size() - 1 ){
			fprintf(rout, ",");
		}
		fprintf(rout, "\n");
	}
	fprintf(rout,"}\n");
};
#endif
/* --------------------------------------------------------------------------------------------------------------*/
bool computestaticpart::mapQ(){
	vector< rank3 > Qmap;
	vector< REAL >  Qvalues;
	vector< rank3 > QImaplocal;
	vector< REAL > QIvalues;
	const unsigned int ksize = kdomain_centered.size();
	for(INDEX ki = 0; ki < static_cast<INDEX>(ksize); ki++){ const REAL k = kdomain_centered[ki];
		int spos = Qvalues.size()-1;
#ifdef DEBUG
		vector< vector<double> > Qmat( numfields, vector<double>( numfields, 0) );
#endif
		vector< vector<double> > Qinv( numfields, vector<double>( numfields, 0 ) );
		for(INDEX ii = 1; ii <=numfields; ii++){ 
			for(INDEX jj = 1; jj <=numfields; jj++){ 
				REAL rval = (* this->func_compute_Q)( *mo, k, rank2(ii,jj) );
				if(rval != 0){
					Qmap.push_back( rank3(ki, ii-1, jj-1) );
					Qvalues.push_back( rval );
				} 
				Qinv[ii-1][jj-1] = rval;
#ifdef DEBUG
				Qmat[ii-1][jj-1] = rval;
#endif
			}}
		// once we have mapped the current matrix for ki
		// we need to invert Q to get QI(k)
		int endpos = Qvalues.size()-1;
		if(endpos - spos > 0){
			// initialize the vector
			if(!invertmatrix( Qinv )){
				return false;
			}
#ifdef DEBUG
			dumpmatrix( Qmat );
			dumpmatrix( Qinv );
			vector< vector<double> > prod;
			multiplymatrix( prod, Qinv, Qmat );
			dumpmatrix( prod );
#endif
			// push_back nonzero values
			for(INDEX ii = 0; ii < numfields; ii++){
				for(INDEX jj = 0; jj < numfields; jj++){
					if(Qinv[ii][jj] != 0.0 ){
						QImaplocal.push_back( rank3( ki, ii, jj ) );
						QIvalues.push_back( Qinv[ii][jj] );
					}
				}
			}
		} 
	}
	this->Qmap->store_coefficients( Qmap, Qvalues );
	this->QImap->store_coefficients( QImaplocal, QIvalues );
	return true;
};
/* --------------------------------------------------------------------------------------------------------------*/
REAL computestaticpart::computeA(const INDEX & ki, const INDEX & pi, const INDEX & qi, const rank3 & ipos) const {
	if( this->func_compute_A ){
		const REAL & k = kdomain_centered[ki];
		const REAL & p = kdomain_centered[pi];
		const REAL & q = kdomain_centered[qi];
		return this->func_compute_A(*mo, k, p, q, ipos); 
	}else{ // computeA using mapped Q, P, and mapped QI
		return computeA_wmap(rank6(ki,pi,qi,ipos[0], ipos[1], ipos[2]), numfields, Qmap, QImap, func_compute_P); 
	}
};
/* --------------------------------------------------------------------------------------------------------------*/
void computestaticpart::mapAhat(REAL * & Afunctions){
	const unsigned int ksize = kdomain_centered.size();
	const int Afunction_size = numfields*numfields*numfields*ksize*ksize*ksize; // dumb guy map, many of these will be unused
	vector< REAL > Ahatvalues;
	vector< rank6 > Amap;
	Afunctions = new REAL[Afunction_size];
	{ 
		REAL fNaN = sqrt(-1);
		for(int x = 0; x < Afunction_size; x++){ Afunctions[x] = fNaN; }
	}
	for(INDEX ki = 0; ki < static_cast<INDEX>(ksize); ki++){ const REAL k = kdomain_centered[ki];
		for(INDEX pi = 0; pi < static_cast<INDEX>(ksize); pi++){ const REAL p = kdomain_centered[pi];
			for(INDEX qi = 0; qi < static_cast<INDEX>(ksize); qi++){ const REAL q = kdomain_centered[qi];
				if( isinsideslot(k,p,q) ) { 
					for(INDEX ii = 1; ii <=numfields; ii++){ 
						for(INDEX jj = 1; jj <=numfields; jj++){ 
							for(INDEX kk = 1; kk <=numfields; kk++){

								rank3 r3(ii,jj,kk);
								REAL Aval = (this->computeA)(ki,pi,qi, r3);
								if(Aval != 0){
									if(findAiter(Afieldindexmap,r3) < 0){
										Afieldindexmap.push_back(r3); // Typically there should be no more than 20 or so of these
									}
								}
								int findex = getAindexfromset( ksize, numfields, ki, pi, qi, ii-1, jj-1, kk-1);
								Afunctions[findex] = Aval;
								assert(findex < Afunction_size);
								if(Aval != 0){
									rank6 r6(ki, pi, qi, ii-1, jj-1, kk-1);
									Ahatvalues.push_back(Aval);
									Amap.push_back(r6);
								}
							}}}
				}
			}}} 
			const_cast<model_class*>(mo)->setAfieldindices(Afieldindexmap);
			Ahatmap->store_coefficients( Amap, Ahatvalues ); 
};
/* --------------------------------------------------------------------------------------------------------------*/
void computestaticpart::mapstaticpart(){
	const model_class & model = *mo;
	const int nksize = this->nkdomain.size();
	// Allocate initial storage for map (overestimate of truly needed size)
	const int prelim_map_size   = nksize * (2*nksize-1) * kdomain_centered.size();
	qgcomplex ** prelim_R_map = new (nothrow) qgcomplex * [prelim_map_size];
	qgcomplex ** prelim_S_map = new (nothrow) qgcomplex * [prelim_map_size];
	rank11 ** prelim_static_map   = new (nothrow) rank11 * [prelim_map_size];
	rassert( prelim_R_map != NULL );
	rassert( prelim_S_map != NULL );
	rassert( prelim_static_map != NULL );
	int * prelim_map_sizes = new (nothrow) int [prelim_map_size];
	rassert( prelim_map_sizes != NULL );
	//  The list of pointers above refers to a list of arrays each of which contains
	//  all the elements for the mk, np, k course grid. 
	// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	// the symmetric part of the integrand depends only on wave-vector magnitude k,p,q, and some terms
	// are re-used, mainly those along the diagonal widths of the slot domain. For purposes of optimization
	// (needed for some domains), the symmetric part is now precomputed in parallel (when available) and stored
	// slots restrictions are used here. To avoid duplicate work, we keep track of the 3 paired indices (ki,pi,qi)
	// done for each course grid point and use it to fill in values if the element was already precomputed. Note that
	// the integrand is fully symmetric in k,p,q, i.e. F(k,p,q) =F(p,q,k) = F(k,q,p) .....
	const unsigned int ksize = kdomain_centered.size();
	const size_t symmetric_data_size = ksize*ksize*ksize;
	REAL * symmetric_data = new (nothrow) REAL[ symmetric_data_size ];
	rank3 * symmetric_data_map = new (nothrow) rank3[symmetric_data_size];
	rassert( symmetric_data != NULL && symmetric_data_map != NULL );
	this->mapsympart(symmetric_data_size, symmetric_data, symmetric_data_map);
	delete [] symmetric_data_map;
	// map Q ....
	if(!mapQ()){
		fatal_stop("Could not invert the Q matrix.");
	}
	// next we map the A map and store it in an array A depends on k,p,q, i,j,k
	REAL * Afunctions = NULL;
	this->mapAhat(Afunctions);

	const int ipairsSize = static_cast<int>(this->Afieldindexmap.size());

	const int cdomainstorage_sizes = ksize*ksize*ipairsSize*ipairsSize;
	// . . . . . . . . . . . . . . . . . . . . . .  
	// Begin construction of static part by looping
	// . . . . . . . . . . . . . . . . . . . . . .
	size_t kiter = 0, citer = 0, piter = 0; // iterate over each k for all sets of mk, np, and k 
	// iterate over all course domain elements (mk, np, k, p, q, i1, j1, k1, i2, j2, k2) 
	for(INDEX mki = 0; mki < static_cast<INDEX>(nkdomain.size()); mki++){ // loop over angular 
		//for(INDEX npi = -static_cast<INDEX>(nkdomain.size()) +1; npi < static_cast<INDEX>(nkdomain.size()); npi++){ 
		for(INDEX npi = -static_cast<INDEX>(nkdomain.size())+1; npi < static_cast<INDEX>(nkdomain.size()); npi++){ 
			// need negative modes too here  
			INDEX np = nkdomain[npi];
			INDEX mk = nkdomain[mki];
			size_t sym_iter = 0;
			for(INDEX ki = 0; ki < static_cast<INDEX>(ksize); ki++){ // loop over each k
				const REAL k = kdomain_centered[ki];
				//const REAL dk = kdomain[ki+1]-kdomain[ki];
				const REAL dk = dkdomain[ki];
				// now declare the domain storage array (size is exaggerated since the domain is actually a slot)
				qgcomplex * cRdomainstorage = new (nothrow) qgcomplex [cdomainstorage_sizes];
				qgcomplex * cSdomainstorage = new (nothrow) qgcomplex [cdomainstorage_sizes];
				rank11 *        c_staticindexmap= new (nothrow) rank11        [cdomainstorage_sizes];
				if( (cRdomainstorage==NULL) || (cSdomainstorage==NULL) || (c_staticindexmap==NULL) ){
					fprintf(stderr, "Problem allocating memory for mki=%i, npi=%i, ki=%i.\n", static_cast<int>(mki),
							static_cast<int>(npi), static_cast<int>(ki) );
					fatal_stop("Memory allocation failed.");
				}
				int cdomain_iter = 0; 
				for(INDEX pi = 0; pi < static_cast<INDEX>(ksize); pi++){ // loop over p
					const REAL p = kdomain_centered[pi];
					//const REAL dp = kdomain[pi+1] - kdomain[pi];
					const REAL dp = dkdomain[pi];
					for(INDEX qi = 0; qi < static_cast<INDEX>(ksize); qi++){ // loop over q
						const REAL q = kdomain_centered[qi];
						const REAL dq = dkdomain[qi];
						const REAL inverse_vol = 1.0/(dp*dq*dk);
						if( isinsideslot(k,p,q) ) {  // only include course elements within the slot
							const REAL & sympart = symmetric_data[sym_iter];
							qgcomplex asympart = asymkpqintegrand( k, p, q, np, mk); 
							// iteration over the field index groups
							for(int iterA1 = 0; iterA1 < ipairsSize; iterA1++){
								for(int iterA2 = 0; iterA2 < ipairsSize; iterA2++){
									const rank3 & A1index = Afieldindexmap[iterA1];
									const rank3 & A2index = Afieldindexmap[iterA2];
									const INDEX & fieldi = A1index[0]-1;
									const INDEX & fieldj = A1index[1]-1;
									const INDEX & fieldk = A1index[2]-1;
									const INDEX & fieldm = A2index[0]-1;
									const INDEX & fielda = A2index[1]-1;
									const INDEX & fieldb = A2index[2]-1;
									int Acommon_index = getAindexfromset( ksize, numfields, ki, pi, qi, fieldi, fieldj, fieldk);
									int A_R_index     = getAindexfromset( ksize, numfields, pi, ki, qi, fieldm, fielda, fieldb);
									int A_S_index     = getAindexfromset( ksize, numfields, ki, pi, qi, fieldm, fielda, fieldb);
									qgcomplex rval = (inverse_vol*sympart*R_coeff*Afunctions[Acommon_index]*Afunctions[A_R_index])*asympart;
									qgcomplex sval = (inverse_vol*sympart*S_coeff*Afunctions[Acommon_index]*Afunctions[A_S_index])*asympart;
									if( (rval != rval) || (sval != sval) ){ // no NaNs allowed
										fprintf(stderr, "Computational sampling occurred outside of allowable domain.\n");
										fprintf(stderr, "Did you turn off the slot domain restriction without using a custom kpqintegrand?\n");
										throw cgcmexception("Fatal error", -1);
									}
									cRdomainstorage[cdomain_iter] = rval;
									cSdomainstorage[cdomain_iter] = sval;
									c_staticindexmap[cdomain_iter]=  rank11( mki, npi, ki, pi, qi, fieldi, fieldj, fieldk, fieldm, fielda, fieldb);
									rassert( cdomain_iter < cdomainstorage_sizes );
									cdomain_iter++;
									citer++;
								}} // close field index loop
							piter++;
						} // close slot domain check
						sym_iter ++;
					}} // close pi, qi loop
					model.printfmsg( "Completed k = %f, kiter = %i, npi=%i, mki=%i\n", k, static_cast<int>(kiter),
							static_cast<int>(npi), static_cast<int>(mki)
						       );
					prelim_static_map[kiter] = c_staticindexmap;
					prelim_map_sizes[kiter] = cdomain_iter;
					prelim_R_map[kiter] = cRdomainstorage;
					prelim_S_map[kiter] = cSdomainstorage;
					kiter++;
			}      // close ki loop
		}} // close angular mode loop

		delete [] Afunctions; // de allocate this memory block
		// write all data back to the R_s static array
		// but first calculate the size
		// model.printfmsg( "Setting Rs sizes.\n");
		int map_size = 0;
		for(int i = 0; i < prelim_map_size; i++){
			map_size += prelim_map_sizes[i];
		}

		model.printfmsg( "Finished setting Rs sizes.\n");
		model.printfmsg( "Allocating Rs array.\n");
		vector< qgcomplex > Rweights( map_size, qgcomplex(0,0) );
		vector< qgcomplex > Sweights( map_size, qgcomplex(0,0) );
		vector< rank11 > static_indexmap( map_size );
		model.printfmsg( "Finished allocating Rs array.\n");
		citer = 0;
		for(int i = 0; i < prelim_map_size; i++){
			for(int j = 0; j < prelim_map_sizes[i]; j++){
				Rweights[citer] = (prelim_R_map[i])[j];
				Sweights[citer] = (prelim_S_map[i])[j]; 
				static_indexmap[citer] = (prelim_static_map[i])[j]; 
				citer++;
			}
			delete [] prelim_R_map[i]; delete [] prelim_S_map[i]; delete [] prelim_static_map[i];
		}  
		delete [] prelim_map_sizes; delete [] prelim_R_map; delete [] prelim_S_map; delete [] symmetric_data;
		model.printfmsg( "Storing static coefficients into container class.\n");
		this->Rstaticweights->store_coefficients( static_indexmap, Rweights );
		this->Sstaticweights->store_coefficients( static_indexmap, Sweights );
		model.printfmsg( "Finished storing static coefficients in container class.\n");
	}
